/*
 * boxmuller.c  Implements the Polar form of the Box-Muller Transformation
 *
 *              (C) Copyright 1994, Everett F. Carter Jr.
 *              Permission is granted by the author to use
 *              this software for any application provided this
 *              copyright notice is preserved.
 */

#include <math.h>
#include <stdio.h>

#include "r250.h"

void boxmuller_init(int seed)
{
    r250_init(seed);
}

double boxmuller(double m, double s)    /* normal random variate generator */
{                                       /* mean m, standard deviation s */
    double        x1, x2, w, y1, v;
    static double y2;
    static int    use_last = 0;

    if (use_last) {
        y1 = y2;
        use_last = 0;
    } else {
        do {
            x1 = (dr250() * 2.0) - 1.0;
            x2 = (dr250() * 2.0) - 1.0;
            w = x1 * x1 + x2 * x2;
        } while (w >= 1.0);
        w = sqrt((-2.0 * log(w)) / w);
        y1 = x1 * w;
        y2 = x2 * w;
        use_last = 1;
    }

    v = (m + y1 * s);

//  printf("boxmuller() -> %f\n", v);

    return v;
}

#ifdef TEST_BOXMULLER

/*
 * Compile using:  gcc -g -DTEST_BOXMULLER -o boxmuller boxmuller.c r250.c randlcg.c -lm
 */
#include <stdio.h>
#include <string.h>

#define TOTAL_BUCKETS   100
#define TOTAL_SAMPLES   100000000

int main(int argc, char **argv)
{
    int i;
    int buckets[TOTAL_BUCKETS];
    unsigned long int value;

    boxmuller_init(1);

    memset(buckets, 0, sizeof(buckets));

    for (i = 0; i < TOTAL_SAMPLES; i++) {
        value = boxmuller(0.5, 0.1) * TOTAL_BUCKETS;

        if (value < 0) {
            printf("Value %d out of range\n", value);
            return -1;
        }

        (buckets[value % TOTAL_BUCKETS])++;
    }

    for (i = 0; i < TOTAL_BUCKETS; i++) {
        printf("bucket[%d] = %d\n", i, buckets[i]);
    }

    return 0;
}

#endif /* TEST_BOXMULLER */
